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The evaluation of Coulomb forces is a difficult task. The summations that are 
involved converge only conditionally and care has to be taken in selecting the appropriate 
procedure to define the limits. The Ewald method is a standard method for obtaining 
Coulomb forces, but this method is rather slow, since it depends on the square of the 
number of atoms in a unit cell. In this paper we have adapted the plane-wise summation 
method for the evaluation of Coulomb forces. The use of this method allows for larger 
computational cells in molecular dynamics calculations. 
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I. Introduction 

When periodic boundary conditions are employed, a molecular dynamics (MD) 
computational cell containing N atoms may be considered to be a unit cell with an 
N-point basis in a crystal lattice. The calculations of long-range Coulomb interactions 
can then be performed with techniques developed for the evaluation of electrostatic en- 
ergies and forces within crystals. While essentially all molecular dynamics simulations 
which include long-range ionic interactions have used the Ewald method 1 to calculate the 
Coulomb energy and forces, alternative methods offer some computational advantages. 
In the following sections, calculations of lattice sums with Ewald, planewise summation, 
and multipole techniques are briefly described and compared. The convergence proper- 
ties of the Ewald summation are shown to be due to effective surface charges included 
implicitly in the sum and a multipole formulation is presented which produces results 
identical to those obtained with the Ewald method. 

II. Lattice summation methods 

In the Ewald method, the potential due to a lattice of unit point positive charges 
is obtained by surrounding each point charge with equal positive and negative Gaussian 
charge distributions and a uniform negative volume charge. For negative charges one 
simply reverses the signs of all charges. The point charges plus the negative Gaussian 
distributions and the positive Gaussian distributions plus the negative volume charges are 
then summed separately, with the convergence of the positive Gaussian-volume charge 
lattice sum improved by performing the sum in reciprocal space. 

When the energy of a collection of positive and negative point charges in the unit 
cell of a neutral crystal is calculated using the Ewald potential, the background charges 
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cancel in the volume of the crystal and the energy sum is generally taken to be that due 
to the point charges alone. 

The overall speed of convergence for Ewald sums can be improved in a simple way 
by adjusting the free parameter in the equations so that the direct and reciprocal sums 
converge at the same rate. Nijboer and de Wette 2 have shown that this will occur if 
a = \/7r/L, where L is the length of a cubic cell edge and a is the free parameter. 
Sangster and Dixon 3 show that the reciprocal sum can be reformulated so that the 
number of terms is proportional to N rather than N 2 , and the free parameter can then 
be adjusted to increase the number of reciprocal space terms and decrease the number of 
direct space terms required to obtain a given degree of convergence. The optimum value 
of a will be that which minimizes the overall computational time. An implementation of 
both the original and modified Ewald methods 4 indicates that the techniques suggested 
by Sangster and Dixon can increase the speed of computation by a factor of about two 
or three. 

A drawback of using the Ewald method for computing Coulomb forces presents 
itself when simulations containing large numbers of particles are considered. For an MD 
cell containing N atoms, the number of Ewald sums required to calculate the forces is pro- 
portional to A 2 . The terms which are summed in an Ewald energy calculation are rather 
complicated, particularly the direct sum. A force calculation requires the calculation and 
summation of terms which are even more complicated, and the repeated calculations of 
lattice sums in a molecular dynamics calculation can be very time consuming, 

An alternative method for performing Coulomb lattice sums is the planewise sum- 
mation method (PSM) first developed by Nijboer and de Wette 5 for dipole lattices and 
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crystals with monoclinic and higher symmetry. In this method, lattice summations are 
computed by applying a two-dimensional Fourier transform over two indices and then 
performing the summation over the third index analytically. There are no additional 
charges introduced to speed the convergence. 

The PSM was later extended to include summation of multipole lattices of all 
orders 6 , summation of multipole lattices in triclinic crystals 7 , and the direct calculation 
of the electrostatic potential 8 . 

The planewise summation method offers some computational advantages over the 
Ewald method. The lattice summation is performed over two indices rather than three, 
there is only a reciprocal space sum and not a direct sum, and the individual terms in 
the sum are generally simpler to calculate than those in an Ewald sum. 

The number of calculations required to find the electrostatic energy by direct 
planewise summation is proportional to N 2 . An additional disadvantage of the direct 
planewise method is the conditional convergence of the lattice sums, leading to a depen- 
dence of the total energy on the choice of the planes in the crystal and inconsistency 
with periodic boundary conditions when the dipole moment of the unit cell is nonzero. 

The conditional convergence of the PSM is due to the general convergence properties 
of the Coulomb lattice sum for a neutral unit cell. The energy sum may be written: 

^2 oo oo oo N N Q'Q' 

E cou i = — H E EI! | z _ ~ 3 _ z | ( x ) 

a=-oo/3=-oc7=-ooi=l j=l I r i ^a(S^ r j \ 

where R a pj = era + (3b + 7c; a, b, and c are Bravais lattice vectors; % is the valence of 
the i th atom; and the prime on the summation over lattice indices indicates omission of 
the terms for which i = j and a = (3 = 7 = 0. 

4 



Using this expression to approximate the energy by evaluation of a finite number of 
terms leads to contradictions with the periodic boundary conditions which are commonly 
imposed on a MD cell. In the usual case of a cell in which the total charge is zero, the 
energy sum will approach a finite limit. However, if the cell has a nonzero dipole moment 
there will be a constant electric field component throughout the MD cell, so that the 
potential at a point on one face of the cell will differ from that at a point translated 
through the cell along a Bravais lattice vector to the opposite face. The existence of a 
non-zero dipole moment also makes the constant electric field component dependent on 
the choice of the unit cell. An example of how this can occur is illustrated in Figure 1. 
Two different, equivalent choices for the unit cell of a tetragonal AB2 ionic compound 
are shown along with the net dipole moment for the unit cell. It is seen that the dipole 
moment changes sign if the locations of the atoms on the corners of the cell are redefined. 
If the derivative of the energy sum is calculated for each of these cells, then each sum will 
contain terms describing an electric field parallel to the dipole moment of the associated 
unit cell, and even though the cells are physically equivalent the electric fields for the 
two cases will differ in sign. 

In an MD simulation with periodic boundary conditions, if a particle leaves the 
unit cell then there is an identical particle which simultaneously enters the cell through 
the opposite face. If the exiting particle is replaced in the energy sum by the entering 
particle, the unit cell is effectively redefined. The configuration of the unit cell could 
change, for example, from one of those shown in Figure 1 to the other in the course of the 
simulation. When this happens, the electric field inside the cell changes discontinuously. 
Since the particles in an MD cell are not taken to have any particular symmetry, a 
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dipole moment is generally present and discontinuities in the electric field and therefore 
the forces is the general rule. 

The dependence of the energy sum on the dipole moment of the unit cell is well 
known. In addition to early studies on the effect of the order of summation of terms 
on the energy totals 9 , the development of fast summation techniques for the Coulomb 
energy of ionic lattices has led to methods which can give different total energy values 
for the same unit cell. The Ewald method and the planewise summation method, for 
example, yield different values if the unit cell has a net dipole moment. 

An alternative to both the Ewald and the planewise method as discussed above 
begins with the separation of the Coulomb potential into multipoles. The resulting 
multipole lattice sums can then be summed by the planewise method. The lattice sums 
do not contain the atomic positions, so if simulations are carried out in a unit cell 
with fixed shape, these sums are constant and need to be calculated only once. For 
simulations in MD cells with variable shape, the lattice sums must be recalculated at 
each time step. In either case, using procedures similar to those of the fast multipole 
method of Greengard and Rokhlin 10 it is possible to construct algorithms for the forces 
which require a number of calculations proportional to N. 

Although a multipole expansion of the Coulomb energy contains the same indeter- 
minacy and dipole dependence of the direct Coulomb sum, a simple modification can be 
made which results in lattice sums which are all absolutely convergent and which gives 
results identical to those obtained by the Ewald method. In the next section, the Ewald 
energy sum is derived as the limit of a finite lattice summation, and the magnitude of the 
implicit surface charges is found. The following section describes the modified multipole 
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formulation and shows the equivalence of this method to the Ewald method. 
III. The Ewald Method for Finite Lattices 

Since Ewald lattice sums are absolutely convergent while direct Coulomb lattice 
sums can be conditionally convergent, there must be contributions to the Ewald sums 
in addition to those from the original point charges and their images. By replacing the 
infinite lattice summation in the Ewald method with a summation over a very large but 
finite lattice, the additional terms can be expressed as surface charges located on the 
faces of the volume occupied by the lattice. 

The Ewald potential includes a sum of terms from point charges, positive and 
negative Gaussian charge distributions, and a uniform background charge. The Gaussian 
distributions cancel at each lattice point, so the electrostatic potential at a point f is 
due only to a lattice of unit point charges plus a neutralizing background charge. This 
potential is: 



M < 1 
$ew(r) = lim V < — 



,. , d 3 r'-^ (2) 
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V is the volume of the unit cell and V^ Q/ g 7 is the volume around the lattice point R a /3j- In 
order to simplify the expression for the potential, the limits of the summations are taken 
to be the same for all three directions. The lattice sum converges with increasing M, 
although the value of the sum depends on the manner in which the limits are approached 
and will be different if we take the limits M a , Mp, andM^ in a different way. 

When the energy of a collection of positive and negative point charges in the unit 
cell of a neutral crystal is calculated using the Ewald potential, the background charges 
cancel in the volume of the crystal. There are, however, regions on the surfaces of a finite 



lattice in which the overlapping volume charges do not cancel. The magnitude of these 
charges is found by consideration of the volume charge contribution to the electrostatic 
energy sum. 

The electrostatic energy in the Ewald formulation is equal to 

g 2 N N 

E ew = tEE $ ^ _J 3') 

i=\ j=l 

= t^coul + ^surf 

e 2 M N N f 1 
Ecoul = -7Z lim Jm'^JTs Z , Z j ( 3 ) 



E w = -- lim y yy M p / dv^ 

The energy due to the uniform volume charge integration can be more simply 
expressed by making a change of variables: 



r = r + ri + 

= eS + V "b + ("c 
n = iia + r]ib + CiC (4) 

= tja + rijb + CjC 

£, T], and ( are the displacements along the a, b, and c directions, respectively. 
With these substitutions, 

, x _p 2 f /-M+i+^i rM+^+iji /-M+i+Ci 1 1 

— / M^oo 2 Y? i J-M-h+ti J-M-\+m J-M-\+Ci l^'-rj-IJ 

(5) 
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With a neutral unit cell, the volume charge is zero within the limits —A to +A, —B 
to +B, and — C to +C. If this volume is subtracted from the sum, the contributions from 
the volume charges are limited to integrations over cells which are at the limits of the 
lattice sums. If the limits are sufficiently large, the volume charge may be approximated 
by a surface charge concentrated at the limit. For example, with M sufficiently large, 



M 



1 1 + Ma- fj | 
-1 



rJk + Ma - fj 



(6) 
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The potential due to the surface charge is found by performing the surface integration 
over rj" and Neglecting edge and corner effects, the total energy due to the surface 
charges is 
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i 

If the surface charge is approximated by a point charge of the appropriate value 
at each surface lattice point, the integrals can be replaced with two-dimensional lattice 
sums. The value of the point charge is found to be numerically equal to that of the 
surface charge. As the limits are taken to infinity, the volume, surface and point charge 
representations become equivalent. 

By inspecting the surface of the finite crystal, it can be seen that the volume charges 
cancel each other exactly only if the net dipole moment of the unit cell is zero. If the 
dipole moment is nonzero, then an Ewald summation contains contributions from surface 
charges with magnitudes given by equation (8) in addition to the direct Coulomb sum 
of equation (1). 

IV. Multipole Lattice Sums 

In terms of the unnormalized spherical harmonics, the inverse distance between two 
points r and f* + R a ^ is 6 : 



V* V* V* V" (- 1 ) l ( k + l - m - n V- r i v * ( ?\ r ik v * ( ^ Y k+l;n+m (R a p 7 ) 

Z^ Z^ Z^ Z^ (U + m\\ ^mV)' 1 kn\' ) D k+l+l 



?-Raf>r-i>\ tor^-kt.rt-1 + + ™)\ ^' ™' 

(9) 



Y lm (f) = P^icosO)^ (10) 

This expression is only valid if r + r' < R a /3y Because of this, it is necessary 
to define a sphere about the origin of the MD cell and restrict the application of the 
multipole expansion to those terms in the lattice sums which come from cells with centers 
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outside the sphere. This radius must be greater than the longest diagonal of the MD cell 
to insure the validity of the multipole expansion in all cases. Contributions from 'near 
neighbors' cells inside the sphere are summed directly, so that the total Coulomb energy 
sum is: 



2 N 

*- = t £ <") 
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+y E E E E E E E 
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With the multipole moments of the MD cell defined as 



N 

Ql m =Y,^rlY lm (r z ) (13) 
i=i 

and the first term in equation (11) defined as E nn , the energy is: 



„ , e 2 ^ ^ ^ ^ Y k+l 

kn Im apy^nn a/3-y 

The lattice sums can be calculated separately for each combination of k + Z and 
n + m. For + / > 4, these sums are absolutely convergent 11 and can be calculated by 
the PSM without ambiguity. The indeterminancy in the total energy is contained in the 
lattice sums with k + I < 4, and since these sums are at best conditionally convergent 
other methods must be used for their evaluation. 

In all cases but one, terms in the sum which contain combinations of k and / with 
k + / < 4 can be shown to be identically zero. If either k or I is equal to zero, then 
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the charge neutrality of the MD cell insures that Qoo = 0. The inversion symmetry of 
the Bravais lattice in combination with the parity of spherical harmonics of odd order 
leads to zero contributions from all odd-valued combinations, including k + I equal to 
1 or 3. The only non-zero terms with k + / < 4 are those with k = 1 and 1 = 1. The 
strength of this term is determined by the collective dipole moment of the charges in the 
MD cell. As expected, if the dipole moment is zero, then Q\ m is zero and there is no 
indeterminancy is the overall sum. 

In order to evaluate the term with k = I = 1 and nonzero dipole moment, a set of 
point charges of arbitrary magnitude are added and subtracted at the center of each face 
of each cell in the lattice. The energy in the MD cell due to these added point charges 
is zero, which in this case is written 

2 M N 
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By combining the contributions from the positive and negative charges at each 
common lattice point, a single lattice sum plus a number of surface terms are created: 



12 



g 2 M N 

2 ^ ' ^ ' - 

a/37=-M i=l 
2 JV M 



i=l 0j=-M 
2 AT M 



e 

' 2> 



2 JV M 

e 



E { 

=a,b,c 










1 r i Ra/3^ — 




2 1 
2 1 


1 ^ " 


~ Ra/3j ~l~ 2 1 


I ^ - 


R-M-iPy + 


a 1 1 
2 1 1 


n - 


Rm/3-y + f 




Qb 






Qb 


I n - 


Ra-M-l-y + 


b | | 
2 1 I 


n - 


RaM~f + 2 1 










Qc 


~| n - 


Raf3-M-1 + 


c 1 
2 1 


1 n - 


Ra0M + § 



-yE* E V ;? r^T r (17) 

^ i=l a /3=-M I r « ~~ KaP-M-l + 5 I I r « ~~ -Ka/3M + 5 I 

The last three terms contain point charges spaced uniformly on the surfaces of the 
finite crystal. Representing these terms by — -E^/ an< ^ performing a multipole expansion 
on the terms in the lattice which are not near neighbors, 
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Referring to the first term in this equation as E nn and the expression in brackets 
in the multipole term as Q% n , this equation becomes: 



- E nn - E^. / + y E E C kn;lmQ*lmQln E ^pfc+Z+l"^ ( 19 ) 
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Values for q a , %, and q c are determined by setting Q\ m = — Qim-The resulting 
values for the compensating charges are 
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Compensating charges for the unit cells of Figure 1 are shown in Figure 2. 

In an orthorhombic lattice the compensating charges are proportional to the dipole 
moments of the unit cell. In triclinic lattices the values are proportional to the projection 
of the total dipole moment along the lattice vectors. Using these values for the com- 
pensating charges, the dipole moment for the MD cell can be equated to an expression 
containing terms which represent a surface charge, a near neighbor sum, and lattice sums 
which are all absolutely convergent: 



2 c ln;lmVlmVln 2^ nk+l+1 
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Here the sums on the RHS are restricted to values of k and I for which k + I > 4. 
When this expression is inserted into the original Coulomb sum, the required form for 
the total electrostatic energy is obtained. 



jp _ TP , E _ E ™P 
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The distance from the original MD cell to the surface charges can be taken to 
be very large compared to the cell dimensions. In this limit the point charges at the 
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summation limits can be well approximated as surface charges with values p a = q a /Ab c , 
Pb = Qb/A ac and p c = q c /A a i,. is the area of the MD cell face defined by lattice 
vectors a and b. A ac and A^ c are defined similarly. 
If the energy of the system is taken as 

E mp = E cou i + E^ r j (23) 

then all of the lattice sums are absolutely convergent and the forces generated by E mp 
are periodic with the lattice. Because of this, the value of the energy is not affected by 
the definition of the unit cell. In fact, E mp is the same as the Ewald energy if the same 
lattice limits are used in both cases, and the surface charge distributions generated by 
each method are identical. 
V. Calculations 

Analytically, the multipole and Ewald methods produce identical numerical results 
when applied to a given MD cell. The choice of which to use in dynamical simulations 
may be made by considering the speed and accuracy of the algorithms available for each 
of these methods. 

The multipole expressions for energy and force were incorporated in a rather simple 
way into energy and force subroutines for molecular dynamics simulations. The near 
neighbor interactions are calculated directly, and the planewise summation method is 
used to calculate the lattice sums for the long-range Coulomb interactions. The resulting 
algorithm includes a number of calculations proportional to N 2 due to the near neighbor 
terms. The multipole calculations also contain a number of terms proportional to N as 
well as a number of calculations proportional to the fourth power of the highest order 
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multipole index included in the sums. 

The multipole expansion is valid for cells with centers more than the maximum 
MD cell diagonal from the the origin, but if this distance is used as a cutoff radius the 
multipole terms will converge very slowly. For maximum efficiency, the number of direct 
calculations within a given cutoff radius must be weighed against the number of terms 
in the multipole expansion which are required to achieve a given accuracy. 

In addition to the error associated with truncation of the multipole expansion, 
additional numerical errors arise from the computation of the lattice sums using the 
planewise summation. Because this method uses Fourier transforms to replace multipole 
lattice sums with more rapidly converging series, it can be applied only to sums over a 
complete lattice. In order to use the planewise method as part of a multipole method, 
the multipole sums are written: 



The primes denotes sums over all lattice points except the origin. The sum over the 
complete lattice can be calculated using the planewise summation method for each value 
of I and m. The terms in the sum over near neighbors are then calculated separately 
and subtracted from the planewise result. 

The subtraction of two lattice sums which are nearly equal in value tends to ex- 
aggerate the numerical errors introduced with the planewise summation method. A 
truncation error in the planewise sums will produce errors relative to the lattice sum 
over all space, which is usually several orders of magnitude larger than a sum which is 
restricted to lattice points outside a cutoff radius. 




(24) 
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For purposes of comparison, subroutines were also written for calculation of the 
Coulomb energy and forces with the Ewald method. The free parameter in the Ewald 
summations was set to a = ^/ir/a, a being the lattice constant. This definition of a is 
used for all calculations described here, with the result that for non-cubic cells the rate 
of convergence depends on the ratios of the lattice constants. 

Error estimates for various sets of computational parameters required to obtain a 
specified error tolerance in the multipole and Ewald methods are listed in Table 1. The 
lattice sums for the planewise and Ewald summations include all terms with the absolute 
value of any summation index less than or equal to Mpsm and M ew , respectively, fi is 
the maximum multipole order included in the multipole sums. 

In order to check the relative efficiency of the two methods, a number of computer 
runs were made with each force routine for the same unit cell. Each run included 81 
force and 3 energy calculations. The parameters used were M ew =3, fj,=14 and Mpsm=7 
for an expected relative error of 1CT 6 . Results are shown in Figure 3. More details are 
found in reference 12. 
V. Discussion 

The multipole subroutine described here includes a number of calculations pro- 
portional to N 2 because of the direct terms in the lattice summation. These terms are 
simpler than those in the Ewald summation, and total run time is less for systems con- 
taining moderately large numbers of atoms. There are many approaches to improving 
the speed of the fast multipole algorithm but these are not pursued here. As imple- 
mented, the multipole routine is faster than the Ewald routine but not overwhelmingly 
so. The relative simplicity of the Ewald formulation makes the Ewald method easier to 
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program and debug. 

In addition to the Coulomb interaction for ionic crystals, it often desirable in simu- 
lations to calculate the energy and forces from interactions which are proportional to any 
negative power of the separation distance. The Ewald method was extended by Nijboer 
and deWette 2 to include this type of lattice sum, and Williams 13 later extended their 
methods to allow multiple atoms in a unit cell. These techniques give formulas which 
can be used to perform fast summation of potentials and forces which would otherwise 
converge very slowly, such as the r -4 dipole-charge interaction commonly encountered in 
ionic models and the r -6 Van der Waals or dispersion interaction. The formulas become 
progressively more complicated as the reciprocal power increases, and it is common to 
use a cutoff radius to calculate these potentials rather than utilize the rapidly converging 
formulation. 

Incorporation of interaction potentials of the form r~ n into a multipole method 
presents some difficulties. The multipole separation used for the Coulomb potential 
cannot be used for the higher power terms; since these terms do not satisfy Laplace's 
equation they cannot be constructed from linear combinations of the solutions to that 
equation. The multipole character of the summations could be preserved by separating 
the lattice sums and particle coordinates through application of a three-dimensional 
Taylor series expansion. These lattice sums could in principle be evaluated with a number 
of operations proportional to N by planewise summation and then differentiated to obtain 
the terms which would be included in the multipole sums. 

An alternative approach for including the higher-power terms in a multipole method 
is found by considering these terms in light of their physical origins. The r -4 term is 
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generally considered to be a charge-dipole interaction, while the r terms usually arise 
from an induced dipole-dipole interaction. If each atom is assigned a polarizability, then 
the dipole moment induced in each atom will be proportional to the electric field at that 
atom. The field at each atom is routinely obtained as part of the force calculation. The 
charge-dipole interaction could then be calculated directly for near neighbors and through 
an additional multipole sum for the long-range contributions. The corrections to the 
electric field at each atom would in turn cause corrections in the atomic dipole moments, 
and repetition of this procedure would be necessary until a specified degree of self- 
consistency is obtained. At each step the only corrections are to the multipole moments of 
the unit cell; the lattice sums remain unchanged throughout the self-consistency routine. 

The direct calculation of dipole interactions through the multipole method would 
have a number of advantages. All of the electrostatic forces due to dipole interactions are 
included automatically; the r -4 and r~ 6 potentials, for instance, need not be considered 
separately. Using fast multipole techniques, it is possible to construct algorithms which 
include all of the electrostatic forces to any required degree of accuracy with a number 
of calculations proportional to N. The short-range forces would be calculated in the last 
step of the fast multipole procedure and would also require a number of calculations 
proportional to N. The method could readily be extended to include quadrupole and 
higher order interactions if required. 
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Table 1. 



Parameters for error tolerances for multipole and Ewald methods. 





n=1.5 

[i Mp Sm 


n=2 

H Mp Sm 


Ewald 

M Ew 


io~ 4 


16 


6 


8 


4 


2 


1(T 5 


20 


8 


12 


6 


3 


io~ 6 


>20 


>8 


14 


7 


3 


10- 7 


>20 


>8 


16 


8 


4 
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Figure Captions 

Figure 1. 

Two equivalent unit cells for tetragonal ZrC>2- Cross-hatched circles are Zr + ions 
and open circles are 2 ~ ions. P is the dipole moment of the unit cell. 

Figure 2. 

Addition of compensating charges which result in zero dipole moments for the unit 
cells of Figure 1. Parallel hatched circles are the compensating charges. 

Figure 3. 

Computer run times for various numbers of atoms per unit cell. Forces are calcu- 
lated by the Ewald method (circles) or by the planewise summation method (triangles). 
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